Resolving singular forces in cavity flow: Multiscale modeling from atoms to 

millimeters 
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A multiscale approach for fluid flow is developed that retains an atomistic description in key 
regions. The method is applied to a classic problem where all scales contribute: The force on a 
moving wall bounding a fluid-filled cavity. Continuum equations predict an infinite force due to 
stress singularities. Following the stress over more than six decades in length in systems with char- 
acteristic scales of millimeters and milliseconds allows us to resolve the singularities and determine 
the force for the first time. The speedup over pure atomistic calculations is more than fourteen 
orders of magnitude. We find a universal dependence on the macroscopic Reynolds number, and 
large atomistic effects that depend on wall velocity and interactions. 
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Processes that span a wide range of length scales pose 
profound theoretical challenges 0, 0- The different 
length scales must be followed with different time res- 
olutions and may require qualitatively different descrip- 
tions of matter. For example, discrete atomistic effects 
may be important in regions of high stress or rapid spa- 
tial variation, while other regions are most naturally and 
efficiently modeled as a continuous medium. Important 
examples of such problems include adhesion and friction 
3, 4], deformation of crystalline solids 0,0, distribution 
and flow of charges at biological interfaces, flow near solid 
surfaces Q, and the many cases where continuum equa- 
tions lead to singularities. 

Several innovative paradigms for bridging between 
atomistic and macroscopic scales have been proposed in 
recent years, and tested against purely atomistic simula- 
tions in small idealized systems B H E El IH III 
Il5j . A few have been applied to specific problems with 
a large range of scales. Notable examples include calcu- 
lations of crack propagation in silicon that include elec- 
tronic structure near the crack tip 0, and calculations 
of indentation with the quasicontinuum method |17| . 
However, these applications have only reached microm- 
eter length and nanosecond time scales, and the main 
effect of large scales is to provide appropriate boundary 
conditions for the atomistic region. 

In this paper we consider a classic problem in fluid 
mechanics where all length scales contribute equally: the 
force on a moving boundary of a fluid filled cavity. Us- 
ing different spatial and temporal resolutions in differ- 
ent regions allows us to treat cavities with dimensions 
of order millimeters and characteristic times of tens of 
milliseconds. Our multiscale approach accelerates the 
calculation by more than fourteen orders of magnitude 
compared to brute force atomistic simulations. Spanning 
many decades in length scale allows us to build a simple 
scaling relation for the total force F that captures both 
atomistic effects and the influence of the only parameter 
in continuum theory, Reynolds number. 

Cavity flow has intrigued scientists because of singu- 



■ 




k, ~:~.z: 






FIG. 1: Geometry and streamlines for cavity flow at Re = 
pUL/n = 6400. a) The square cavity has edge L = 10 6 er ~ 
0.3mm in the x — y plane, the top wall moves to the right at 
speed U, the flow is independent of z, and the natural time 
scale to reach steady state flow is ~ 40ms. Solid (dash-dotted) 
streamlines indicate clockwise (counter-clockwise) flow. In 
(b-d), flow near each corner is magnified by the indicated 
factor. The coupling of different resolutions is illustrated in b. 
The coarser solution provides boundary conditions along the 
outer boundary of the finer grid (dashed lines) and obtains 
them along the dotted lines. In (d), dotted lines indicate 
the boundaries of the ~ 12a wide regions that are treated 
atomistically. 



larities that can not be resolved by purely continuum 
methods. Figure QJi illustrates the cavity geometry. The 
top wall is displaced to the right at fixed velocity U and 
the other walls are stationary. The traditional contin- 
uum approach to this problem uses the Navier-Stokes 
(NS) equations and no-slip boundary conditions at the 
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walls [l8L [Tgl l2*cj . The no-slip condition requires that the 
fluid velocity u vanishes near static walls, and equals U 
at the moving wall. The discontinuity in boundary con- 
dition at corners between moving and fixed walls causes 
the stress to diverge as the inverse of the distance from 
the corner, r. The total force on the wall is the integral 
of the stress, and diverges logarithmically in continuum 
theory. 

Koplik and Banavar pioneered the use of molecular 
dynamics (MD) to study the singularities in cavity flow 
[211 ] . They observed a breakdown of the no-slip boundary 
condition within atomic distances from the corner, and 
a corresponding saturation of the stress. Similar effects 
cut off stress singularities in the closely related problem 
of spreading fluids j2^,|2^|. Koplik and Banavar's purely 
atomistic approach limited the cavity length to L = 17 a, 
where a ~ 0.3nm corresponds to a molecular diameter. 
We recently extended L to 250c ~ 75nm using a hybrid 
method that treated singular regions atomistically and 
the remainder of the cavity as a continuum. This ap- 
proach allowed us to analyze the breakdown of the con- 
tinuum boundary conditions over a wider range of U , and 
to determine its microscopic origins . The stress devi- 
ates from the singular continuum solution for r < S(U). 
At low U, the intrinsic discreteness of atomic fluids leads 
to S ~ <7. At large U, the interfacial stress is high enough 
to produce non-Newtonian effects and S rises linearly 
with U. 

To study cavities with millimeter scale dimensions, two 
major improvements must be made on previous work. 
The first is to vary the resolution in the continuum re- 
gion to efficiently describe the rapid increase in velocity 
gradients near the singular corners. We chose to do this 
using a local refinement approach (2^. The second is to 
span the wide range of time scales associated with dif- 
ferent spatial resolutions. While the motion of atoms 
near the corner must be followed with time steps of or- 
der 10 _14 s, the time for flows to equilibrate on millimeter 
scales, ~ 100L/U, may be of order seconds. To overcome 
this obstacle we use the optimum time step to obtain the 
steady state flow at each spatial scale, and then enforce 
self-consistency between scales. 

The details of our approach and a sample flow field 
are illustrated in Fig. ^ At the continuum scale the flow 
is independent of z, and satisfies the NS equations with 
vicosity ji, fixed density p, and no-slip boundary condi- 
tions. At each scale the NS equations are discretized on 
a square grid of cells with width h, and the steady state 
flow is obtained using an artificial compressibility method 
[26| with time step of 1/4 - 1/2 h/U. On the coarsest 
scale h = L/256. This resolution is inadequate near the 
corners, where h is decreased by successive factors of two 
using an iterative refinement scheme. 

At each stage of the iteration, solutions at two reso- 
lutions provide boundary conditions for each other. The 
geometry is illustrated in Fig. QJ). Both resolutions use a 
64 by 64 array of square cells. The finer grid lies in the in- 
ner quarter of the coarse grid (dashed lines) , and receives 



boundary conditions on its outer edge. It in turn provides 
boundary conditions for the coarser grid along the dot- 
ted lines. The overlap region between dashed and dotted 
lines prevents discontinuities due to sudden changes in 
resolution |9L ITol[l2| . 

This refinement scheme is iterated until the overlap 
region reaches nanometer scales. There the finest res- 
olution results are obtained from MD simulations that 
can be extended all the way into the corner (Fig. QJi). 
At the outer edge of the MD region the mean atomic 
velocity is constrained to follow the finest continuum so- 
lution (h = 0.95tr) and particles are added or removed to 
match the continuum flux dSH. Average MD velocities 
provide boundary conditions for the continuum solution 
along an inner square whose edge is 6 cells long. A global 
steady state solution is obtained by iterating from coars- 
est to finest scale and then back to the coarsest scale until 
all boundary conditions are consistent. This typically re- 
quires ten to twenty iterations, depending on the desired 
accuracy. 

To obtain a smooth solution, the continuum model 
must accurately describe the atomistic behavior at the 
outer boundary of the MD region. This requires consis- 
tent choices of n and molecular interactions. Following 
previous work |2lLl24| . we consider fluid atoms of mass m 
interacting with a Lennard-Jones (L J) potential of char- 
acteristic energy e and diameter a. The potential is trun- 
cated at r c = 2.2cr, and the mass density p = 0.81to(j~ 3 . 
The geometry and interactions of the crystalline walls are 
chosen to produce a no-slip boundary condition far from 
the corner [24|. Discrete wall atoms are on the sites of 
a (111) surface of an fee crystal of lattice constant 1.204 
a and interact with the fluid with a LJ potential with 
energy e wf = 0.95e. 

Within the MD region, the motion of particles is fully 
three-dimensional. However the mean flow velocities are 
independent of z, and periodic boundary conditions with 
period L z are applied in this direction. The equations 
of motion are integrated using the Verlet scheme with 
time step AtMD = 0.005tLj, where t^j = (mcr 2 /e) 1 / 2 is 
the characteristic time of the LJ potential. A constant 
temperature ksT = l.le is maintained by applying a 
Langevin thermostat |27| with dampin g ra te T = UJ^j in 
the z direction. The dynamic viscosity 10J of the LJ fluid 
fi = 2.14eiijcr -3 is used in the NS equations. 

Two lengths characterize transitions in flow behav- 
ior near each corner. Incrtial effects are significant for 
r > rj = fi/ pU i while deviations from continuum behav- 
ior occur for r < S(U). At intermediate scales our results 
follow the analytic solution for Stokes (non-inertial) flow 
[l8j | . Near the top corners the streamlines are scale invari- 
ant because u only depends on the angle relative to the 
moving wall. The stress diverges as 1/r since u changes 
by U over a length of order r. A series of counter-rotating 
vortices forms near the bottom corners. The change in 
scale between panels in Fig. Q]was chosen to illustrate the 
predicted self-similarity under r — > r/16.4 and a change 
in direction of rotation 0, [2^ . The vortices are cut off 
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FIG. 2: Shear stress as a function of distance r from the 
upper left corner for U = 0.27a /tij . Qualitatively similar 
scaling occurs near the upper-right corner. In (a) Re = 6400, 
L = 6.25 x 10 4 cr, ri = 9.8a, S fa 2.2a, lines of alternating 
color indicate continuum results from successive resolutions, 
and asterisks show MD results. Results for Re = 25 to 6400 
obtained by increasing L by factors of two are compared in 
(b). 
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FIG. 3: Dimensionless force per unit length on the sliding 
wall F/pU as a function of Reynolds number Re — pUL/p. 
Results for U = 0.27a/t LJ (circles) and U = 0.06Sa/t LJ 
(squares) follow Eq. (dashed lines) with S — 2.2 and 0.55 
respectively. The inset compares the calculated force vs. ve- 
locity at Re = 400 (circles) to Eq. [I] (line). 



at S(0) ~ a by atomistic effects, and this is responsible 
for the small deviation between panels c and d. 

Fig. |2I a) shows a plot of the shear stress t along the 
moving wall at Reynolds number Re = pUL/p = 6400 
as a function of distance r from the upper-left corner. 
Alternating colors are used for data from different res- 
olutions to illustrate the smooth matching. The as- 
terisks are from the force per unit area on wall atoms 
in the atomistic region. At intermediate scales, the 
stress follows the divergence predicted by the Stokes so- 
lution r = (47r)/(7r 2 - 4)(jiU)/r (straight line). At large 
scales the stress decays more slowly than 1/r, and at 
r < S w 2.2ct the stress singularity is cutoff by atomistic 
effects. Fig. [21 b) shows that changing Re by increasing 
L only changes the flow in the outer region. At scales of 
order L/2 the other corners become important, but for 
r <C L/2 all results fall onto a common curve. This shows 
that S is independent of L and only depends on U and 
atomic properties. 

The total shear force on the moving wall is an integral 
of the shear stress along it. Fig. [2] shows the force per 
unit length along the z— direction, F, as a function of 
Re. If the Stokes solution applied, each factor of two in 
length scale would contribute the same force and the total 
would diverge. This divergence is cutoff by atomistic 
effects at r < S but inertial effects enhance the force from 
r > rj. To determine F, one must resolve the inertial 
effects between r% and L = Re n as well as the atomistic 
behavior at the corner. Our approach enables us to span 
this wide range of scales (> 10 5 ) for the first time. 

If atomistic effects were not important, the dimension- 
less force on the wall / = F / pU would only depend on 



Re. However, Fig. [21 shows that changing U at fixed 
Re produces large shifts in /. A complete description 
of / can be obtained by considering the separate con- 
tributions, fi, of the three scaling regimes in the stress 
shown in Fig. [21 They are well-separated as long as 
Re = L/rj 3> 1 and the Reynolds number at molecular 
scales Re m = S/rj <C 1. 

The dependence on Re comes entirely from the outer 
region r > 77 , which is why curves for different U in Fig. 
[21 are related by a constant shift. As shown in Fig. [2[b), 
increasing Re extends the range of the inertial tail that 
contributes to the dimensionless force, causing fs(Re) to 
rise with Re. For Re > 10 our numerical results are well 
described by f 3 = a + bRe a with a = 3.85, b = 1.98 and 
a = 0.434. 

In the inner region near each corner, r < S, the stress 
is determined entirely by U and atomic properties, such 
as tLj, a, p, and wall geometry (Fig. Gib)). For fixed in- 
teractions, we can write the contribution from this regime 
as fi(UtLj/a). In fact we find /1 has a constant value 
of about 4.3. The reason is that the range of integration 
grows linearly with S while the stress scales inversely 
with S. Indeed assuming the stress for r < S is equal to 
the Stokes stress at S yields h = 8n/(7r 2 - 4) w 4.28, 
which is consistent with the numerical results. The value 
of /1 is as much as 20% of the total wall force in Fig. [21 
yet it comes from an inner region that is never larger than 
a few molecular diameters. This is clear evidence of the 
direct effect of atomic scale behavior on the macroscopic 
force. 

In the intermediate region, S < r < r\r, the Stokes 
solution applies, and can be integrated to give f%{Ti / S) = 
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8ir/(ir 2 — 4)ln(rj/S). This term is responsible for the 
velocity dependence in the total dimcnsionless force: 

4j = fi + -r L -; l ^ri/S) + a + bRe a . (1) 

The ratio S/rj reflects the range at which atomistic ef- 
fects cut off the Stokes region, and decreases as U de- 
creases. The inset to Fig. shows that numerical re- 
sults for the velocity dependence of F are consistent with 
S = Sq + kUt^j with So = 0.3cr and k = 7. This form 
for S is consistent with our previous studies |24| . which 
showed S approached a constant of atomic scale at low 
velocities and rose linearly with U at high velocities due 
to non-Newtonian effects. 

Our multiscale method has allowed us to span the wide 
range of length scales (> 10 5 ) that contribute to the drag 
force on the moving wall of an ideal cavity, and to extract 
a simple and accurate physical description (Eq. ^) of the 
important contributions from each scale. We have used 
it to treat cavities with dimensions of order millimeters 



and natural time scales approaching seconds, and the ap- 
proach is readily extended to still larger scales. Its ma- 
jor limitation is that it assumes a steady state solution, 
while cavity flow becomes turbulent at Re > 8000 [2^| . 
It may be possible to extend simulations into the turbu- 
lent regime using ideas like the "equation-free" approach 
of Kevrekidis et al.^5|. Turbulent fluctuations occur on 
the longer time scales associated with coarser resolutions, 
and the finer scales can be iterated to a quasistatic state 
that follows the coarse solution. This would allow cal- 
culations with the same coarse time step that is used in 
a completely continuum description, while retaining the 
crucial atomic detail. We hope that our work encour- 
ages such efforts, as well as applications to other impor- 
tant multiscale problems such as contact-line motion and 
contact mechanics. 
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